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Abstract 

We report Brownian dynamics simulation of the out-of-equilibrium dynamics (aging) in a col- 
loidal suspension composed of rigid charged disks, one possible model for Laponite, a synthetic clay 
deeply investigated in the last few years by means of various experimental techniques. At variance 
with previous numerical investigations, mainly focusing on static structure and equilibrium dynam- 
ics, we explore the out-of-equilibrium aging dynamics. We analyze the wave-vector and waiting 
time dependence of the dynamics, focusing on the single-particle and collective density fluctuations 
(intermediate scattering functions), the mean squared displacement and the rotational dynamics. 
Our findings confirm the complexity of the out-of-equilibrium dynamical behavior of this class of 
colloidal suspensions and suggest that an arrested disordered state driven by a repulsive Yukawa 
potential, i.e., a Wigner glass can be observed in this model. 
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I. INTRODUCTION 



The out-of-equilibrium dynamical behavior of soft matter — and gels in particular — still 
constitutes an interesting puzzle in condensed matter physics [HQ, sj]- Indeed, despite the 
fact that out-of-equilibrium non-ergodic soft-materials are ubiquitous in the everyday life, 
the mechanisms governing their dynamics are still unclear. An adequate comprehension of 
the above mechanisms is important for at least two main general reasons. First, nano-sized 
particles ranging from hundreds of nanometers to microns are currently used as building 
blocks of nanostructered composite materials, therefore having strong implications in indus- 
trial processing and technology development. Second, soft matter systems can help with 
fundamental questions about the nature of condensed matter. Indeed, due the fact that 
interactions can be chemically tuned (strength, range) with a high degree of accuracy, they 
can provide a good implementation of particular limits as, for instance, the hard sphere 
purely repulsive potential, and used to test most part of the theoretical models developed 
to address the physics of liquids. 

From an experimental point of view, relevant time and length scales in colloidal systems 
are more easily accessible than in the case of simple atomic and molecular supercooled liq- 
uids and glasses Q. Nowadays dynamics can be studied to relatively long time scales and 
average collective dynamics is accessible together with single-particle behavior. Different 
spectroscopies are commonly used, ranging from light and X-ray scattering to confocal mi- 
croscopy Q,y,lf3]- These techniques allow one to access the instantaneous state of the system 
not only under the form of spectra or correlation functions of the observables of interest, 
but also as direct snapshots (configurations) of the system. The direct knowledge of the 
space position of the particles can be directly used to calculate all the relevant quantities or 
also for visual inspection [7|, sometimes giving a direct physical intuition of the mechanisms 
controlling the dynamics. 

While technological improvement has strongly developed the above experimental scat- 
tering techniques, theoretical descriptions often lay behind, due to difficulties in modeling 
the interparticle potentials and in a proper description of the role of the solvent. Indeed, 
properties of colloidal systems have often been studied by classical Newtonian molecular dy- 
namics simulations, where the interaction among particles (either spherical or anisotropic) 
are schematized by simple isotropic potentials. The above approach completely neglects 
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the interaction of the colloidal particles with the solvent where they are dispersed, not to 
mention hydrodynamic effects. Important steps have been made in the direction of more 
realistic descriptions implementing Lattice Boltzmann and Fluid Particle Dynamics meth- 
ods where hydrodynamic effects are also considered ^, 0, Q| • These approaches are very 
demanding from a computational point of view, strongly limiting both the size of the consid- 
ered systems and the time and length scales accessible by the simulation. An intermediate 
approach is provided by Brownian dynamics simulation (i.e., neglecting hydrodynamic in- 
teractions) by correctly taking into account the interaction of a possibly anisotropic colloidal 
particle with the solvent and in the presence of a more realistic interaction potential among 
the constituents. This intermediate approach, being computationally feasible, could help in 
understanding out-of-equilibrium dynamics on very long time scales. 

In this work we introduce a new Brownian dynamics algorithm and study the out-of- 
equilibrium dynamics of a colloidal suspension formed by strongly anisotropic discotic-like 
molecules, which is believed to be a good model for Laponite. Laponite is an interesting 
material; it is a synthetic cl ay, deeply inve stig ated in the last few years by means of several 



experimental techniques 



c cl ay deeply investigate 

HHHQ,Q,yii2 



171 ] . This system presents many different o pen 



problems, ranging from the correct description of the phase diagram 



to the out-of-equilibrium and aging dynamics |23j, |2J, |25(, also under particular external 



ts many aitiere 
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conditions |26j|. Previous computer simulations have focused on the static structure and 
the role played by the distribution of the charges on the platelet HHQ- A Brownian 
dynamics simulation has been presented in Ref. |3J|, with a focus on structure and dynamics 
at equilibrium. Here we present an analysis of the out-of-equilibrium dynamics of a model 
introduced in Ref. j^], where a Laponite platelet is modeled by identical negatively charged 
interacting sites, uniformly distributed in a two-dimensional disk-like geometry. 

The paper is organized as follows. Section |TI] contains an overview of the experimental 
phase diagram of Laponite and a description of experiments pointing to the existence of an 
ergodic to non-ergodic phase transition. Section lTTTl is a description of the model we have used 
to describe the interaction between two Laponite platelets, while in Section IIVI we briefly 
recall the principles of Brownian dynamics and introduce our algorithm for the integration 
of the equations of motion in the case of strongly anisotropic colloidal molecules. Section IS 
reports the details of the computer simulations. The next four Sections contain a detailed 
investigation of the long-time aging dynamics of the model. Section IVII reports a study of 
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the mean-squared displacement of the molecules; the one-particle and collective dynamics of 
density fluctuations are included in Sections IVIII and rVHIl respectively. Section llXl concludes 
the study of the out-of-equilibrium dynamics with the rotational dynamics behavior of the 
system. Finally, Section |X] contains a discussion of our main results and the conclusions. 



II. LAPONITE 



Laponite is a very interesting material. Not only it is used in several industrial pro- 
cesses, but it has also been widely investigated by means of the most powerful experimen- 
tal techniques. It is a synthetic clay, formed by very thin cylindrical platelets of radius 
r a ~ 12.5 nm and thickness d ~ 1 nm. The total uniform surface charge is Q ~ Ze, with 



Z ~ —700. Laponite phase diagram has been investigated in details in Ref. 18|. Recently 
significant modifications to the phase diagram of Ref. Hi have been proposed by different 
authors HQ 3 2l|- 

In a series of papers H? 12, Q, 14], dynamic arrest in Laponite has been interpreted 



as formation of a Wigner glass. In particular, aging of a glass is related to the presence of 



cages of particles at high concentration, mainly due to repulsive interactions, w 



rile gelation 



corresponds to cluster formation due to attraction. Recently Ruzicka et al. |19l . l20l |2l( have 
performed dynamic light scattering experiments to study aging and structural arrest for both 
low and high concentrations. Measurements performed over long periods of time have shown 
two different routes to arrest, depending on the role of attractive short-ranged interactions. 



The situation is reminiscent of t 
gives rise to Wigner-type glasses 



le case where competing interactions are present, which 
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171 ] based 



A different interpretation has been proposed by Nicolai and co-workers 
on a picture in which aggregation between the platelets is induced by salt addition, which 
reduces electrostatic repulsion and produces fractal aggregates that will eventually form a 
gel. It was suggested that attractive interaction leads to aggregation of the Laponite particles 
and the formation of a house of cards structure. It was recently shown that the presence 
of positive c harg es on the rim of the Laponite disks is necessary to induce aggregation and 



gelation 



ve cnarg 

HQ- 



These charges were neutralized by added pyrophosphate and aggregation 



and gelation was slowed down, even though the resulting ionic strength was increased. 
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III. THE ATOMISTIC MODEL 

We consider here the realistic interaction model introduced in Ref. |27[ and conventionally 
named "model A ", where the total negative charge is uniformly distributed over the surface 
area and the possible presence of positive rim charges is neglected. In this respect, this 
work constitutes an effort in the direction of quantifying the routes to dynamic arrest driven 
by repulsive interactions which progressively establishes via structural reorganization of the 
system. Future work must address the issue of the differences in arrest introduced by the 



presence of attractions mediated 
Following the authors of Ref. 



3V 



oppositely charged rim particles 



271 ] . each Laponite platelet is schematized as a rigid disk, 
formed by v = 61 charge sites disposed on a regular mesh with grid points spaced by r D /4. 
Each site carries a charge q = Qjv ~ — 11.47e. The solvent (water in the present case) is 
treated as a continuum of dielectric constant e = 78. Three among the v charge sites are the 
dynamical sites, i.e., the ones whose dynamics is explicitely integrated in time to generate 
the Brownian molecular dynamics trajectory, and have mass m = 747900 m p (here m v is the 
mass of the proton). The remaining interaction sites are mass-less and they are only taken 
into account in the calculation of the forces acting on the platelets. Their positions can be 
trivially calculated at any time in terms of the coordinates of the dynamical sites. 

The total interaction energy between two platelets is the sum of v 2 site-site screened 
Coulomb electrostatic interactions of the Yukawa form. The interaction between sites k and 
I at positions fk and r\ and pertaining to two different platelets is therefore of the DLVO 



form 



where r^i = ||rfe — and Xd is the Debye screening length of the micro ions, i. e., 



Xd y 47r(n + 4 + n_zi)e 2 ' (2) 

In the above equation, (n + ,z + ) and (n_,z_) are the concentrations and the valences of 
positive counterions and negative coions coming from the platelets and the salt added to the 
solution. 
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IV. BROWNIAN DYNAMICS AND THE INTEGRATION ALGORITHM 



In this article we introduce a novel algorithm to integrate the Brownian equations of 
motion which, although neglecting hydrodynamics effects, takes into account the presence of 
the solvent in the dynamics. Most Brownian dynamics simulations are based on the original 



Ermak's algorithm 36] developed for atomic systems, where the friction coefficient associated 



to the damping force and acting on spherical objects is isotropic. The Ermak's algorithm 



371 ] to the case of molecules formed 



has been generalized by van Gunsteren and Berendsen 
by spherical centers of interaction, for which the choice of an isotropic friction coefficient is 
still appropriate. In contrast, in the present case of a strongly anisotropic rigid molecule, 
one should take into account the fact that the platelet moves differently in the directions 
perpendicular and parallel to its symmetry axis. 

To the best of our knowledge, the only other attempt to take into account in a reasonable 
fashion the anisotropy of a discotic-like particle in a Brownian dynamics computer simulation 
has been the one by Odriozola et al. , based on the integration of the Langevin equations 
of a rigid body, separately for the center of mass and for the orientation. In the present 
work we prefer to rely on the original work by Ermak, which allows us to find an optimal 
balance between a reasonable realism of the molecular model and an acceptable efficiency. 
This gives us the possibility to follow the dynamics of the system on time scales much longer 
than the ones investigated so far. 

In our algorithm the viscous damping and the Brownian forces act on each of the three 
sites with mass but the friction coefficient, £, is chosen different in the directions perpendic- 
ular and parallel to the symmetry axis of the platelet. The numerical values of the 
two friction coefficients must be chosen properly, in order to reproduce a realistic dynamics 
for the single platelet, as we will describe in details in Section El £j_ and £y are used to 
evaluate the following numerical coefficients: 

1 — r a 1 — r a 

dg = e<** c« = i-^ ^ = V^ (3a) 

„. 2 _ ksT 6_t_ ( _ 3 -4e-^ + e- 2 ^ _ k B T _ 5 



w =^r a { 2 ut — ) K) = ^r (1 - e ] (3b) 



k B T 1_ 



;i - e-*- at ) 2 , (3c) 
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where a =_L, || and St is the integration time step. 

For each dynamical site of mass m of the platelet, and at the time (t + St), the integration 
algorithm can be schematized as follows: 

1. Each pair of components of the vectors 

Sr = (Sr x ,Sr y ,Sr z ) (4a) 



Sv = (Sv x ,Sv y ,Sv z ), (4b) 

is sampled from a bivariate Gaussian distribution with zero mean values, variances 
given by (er" ) 2 and (v?) 2 (Eq. (|3b"jl) and a correlation coefficient determined by c! 



(Eq. (jSc))) (see Ref. 38]). With this step we produce a realization of the stochastic 



process appropriate to generate Brownian motion 

2. Transform position, r(t), and velocity, v(t) evaluated in the laboratory-fixed frame 
at time t, to r"(t) = (r' x ,r' y ,r z ) and v'(t) = (v' x ,v y ,v' z ) in the body (platelet)-fixed 
reference frame, where the z-axis is chosen parallel to the platelet symmetry axis 

r f (t)=Rr(t) (5a) 



v , (t) = B,v(t); (5b) 

here R is the appropriate orthogonal matrix corresponding to the change of the refer- 
ence frame 

3. Update positions in the body-fixed reference frame to (t+St) and velocities to (t+St/2), 
according to Ermak's prescription: 

r' (t + St) = r p '(t) + (%v' p 8t + c^^&St 2 + Sr p (6a) 

v' p (t + |) = c«v'p(t) + (c? - c«)^St + Sv,, (6b) 

where a =\\ if (3 = x,y, and a =± if /3 = z. Note that Fp(t) are the forces evaluated 
at time t and also include the intra-platelet contributions coming from the mass-less 
sites 
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4. Transform positions and velocities to the laboratory-fixed reference frame: 

r%t + 6t)=R- 1 r f (t + 6t) (7a) 



v(t + St) = R" V(t + St) (7b) 

I — I 

5. Adjust positions and velocities using the RATTLE algorithm [38j, in order to fulfill 
the rigid constraints (fixed distances) among the three sites with mass. 

6. Evaluate the forces Fp(t + St) acting on the site at time (t + St). 

7. Update velocities to t + St following Ermak's scheme: 

vp(t + St) = v p {t + ^) + <%5t Fp{t + St \ (8) 
where, again, a —\\, if /5 = x, y and a =_L, if (3 = z. 

8. Finally, adjust velocities using the RATTLE algorithm to assure zero relative velocity 
among the dynamical sites. 

V. SIMULATION DETAILS 

We have simulated a system composed of N = 108 Laponite platelets of radius r = 
12.5 nm. The time unit is to = 73.6 ps and the time step used to integrate the equation of 
motions is St = 0.1 to- The system is initially thermalized at a high temperature, T = 500 K, 
and at time t w = 0, it is instantaneously quenched at the room temperature, T = 300 K. t w is 
the waiting time, i.e., the elapsed time after the quench, and it is an additional relevant time 
scale for the dynamics of the system, as we will see below. The total volume of the system 
is V = 1.515 x 10 6 nm 3 , which corresponds to a volume fraction <p = 3.5% in the case we 
assume that each platelet occupies a volume equal to 47rr^<i = 1963.5 nm 3 . The side length 
of the simulation box is L = 114.84 nm and cubic periodic boundary conditions are used. 
The Debye screening length value chosen for the interaction potential Eq. ((TJ is A n = 3 nm, 



a value which is believed to be appropriate to realistic experimental conditions [2jj. The 
interaction potential is truncated when the distance between the two interacting sites exceeds 
a cut-off radius r r = 10 nm. 
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The numerical values of the friction coefficients, £j_ and £y, defined above have been chosen 
according to the following argument. A Laponite platelet can be assimilated to an oblate 
ellipsoid of semi-axes (a,b,c), with a = b and a > c, which diffuses in a solvent (water) 
characterized by a viscosity rj. In the case of an oblate ellipsoid, the translational diffusion 
coefficients, D\ and DT — corresponding respectively to displacements perpendicular and 
parallel to the platelet surface — , and the corresponding rotational diffusion coefficients, 
and — that correspond to rotations about axes perpendicular and parallel to the 
platelet surface respectively-, ean be evaluated exactly (see, for iustauce, Ref. Q). More 
specifically, assuming that: i) we consider as a solvent pure water, i.e., r\ = 1.002 cP; ii) the 
freely diffusing oblate ellipsoid has an aspect ratio similar to that of the Laponite platelets, 
i.e., a = 12.5 nm and b = 0.5 nm; Hi) the temperature is T = 300 K; we have chosen the 
friction coefficients £j_ = 0.016 ps _1 and £|| = 0.063 ps~ x . 

As we will see below, following the temperature jump, the simulated system does not 
reach equilibrium on the time scale of the simulation. This implies that two-time correlation 
functions will depend separately on the time difference t and the time t w elapsed from the 
quench. Hence, average over starting times must be substituted with an average over an 
ensemble of several statistically independent initial configurations. We have considered 100 
independent configurations of the systems which have been equilibrated at high temperature 
T = 500 K and, at time t w = 0, quenched instantaneously to ambient temperature T = 
300 K. We have calculated the Brownian dynamics trajectories for each sample for 10 6 
integration time steps, corresponding to a total time interval of 7/is. The total CPU time 
amounts to about 14000 hours on a farm of 10 64MHz Athlon CPUs. The evolution of 
the 100 independent systems has been recorded and system configurations have been stored 
logarithmically spaced in time for the analysis presented in the following Sections. 

Following the temperature jump from high temperature, the energy of the system relaxes 
but never reaches thermodynamic equilibrium at the lower temperature. In Fig. ^ we show 
the energy relaxation for each sample (solid lines) and the average value (open circles). Note 
that, due to the repulsive nature of the interaction potential, the total potential energy is 
always positive. The relaxation follows a power law of exponent 7 ~ 0.35 and the equilibrium 
value (estimated at ~ 206.6 eV) is not reached on the simulation time scale. Hence, the 
system is always out-of-equilibrium, and explores a series of metastable states of lower and 
lower energy. In what follows we will characterize the aging dynamics, and show that the 
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slowing down of the dynamics strongly depends on the waiting time, t w . 

In the main panel of Fig. El we show the static structure factor, S(Q) = (p(Q)p*(Q)), 
where p(Q) = ^/M^f 1 exp(i Q-X{) is the density fluctuation of wave-vector Q. Here and in 
the following, the coordinates Xi indicate either the center of mass of platelet Ri (in this case 
M = N) or the position r*j of all sites (in this case M = 61 x N). The figure shows our results 
for both centers of mass (open symbols) and sites (closed symbols). In particular, latter data 
contain detailed informations on the form factor of the clay platelets and, therefore, they 
are the results most appropriate for a comparison with experimental measurements. We 
note, however, that the smallest wave-vector value accessible in our simulation has modulus 
Q m = 2n / L ~ 0.055 nm _1 , obviously far from the values typical of the most popular 
scattering techniques. Each point shown is a spherical average in momentum space with a 
resolution AQ = 0.05 nm" 1 and up to 200 vectors have been considered for each value of 
Q. The data have been calculated from the configurations produced at times larger than 
10 6 ps and no changes have been found in the static structure at later times. Only at much 
shorter times small changes are evident in the structure of the second sub-peak for the case 
of interaction sites, as shown in the inset of Fig. El 

In summary, from the above data it is evident that, following the temperature jump, 
the system starts to explore out-of-equilibrium states of lower and lower energy. During 
the above process the static structure does not significantly change at long times, as also 
found in the case of molecular glass formers. This is at strong variance with the long-time 
structural dynamics — as shown by the mean-squared displacement and density fluctuations 
correlation functions — , which show spectacular changes, as we will report in the following 
Sections. 

VI. THE MEAN-SQUARED DISPLACEMENT 

In this Section we analyze the centers of mass mean-squared displacement of the center 
of mass of the platelets, 

1 N 

(r 2 (t, t w )) = - J2(\ R k(t + tw) - R k {t w )\ 2 )- (9) 

k=l 

here () is the average over the initial conditions. Note that in the calculations we have 
taken care of subtracting the displacement of the center of mass of the whole system, which 
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does not vanish in Brownian simulations, to reveal the intrinsic dynamics of the particles. 
Indeed, in the case where the dynamics is very slow (as it is the case here), at long times 
the diffusion motion of the total center of mass can be relevant. In Fig. |H]we plot the mean- 
squared displacement at the indicated values of t w . Following an initial transient regime, the 
curves reach a plateau, reminiscent of the cage effect in glassy systems, and whose lifetime 
increases on increasing t w ; only at later times, eventually, the curves move from the plateau 
and apparently bend toward a constant value. From these data we can tentatively argue 
that, for values of t w longer than the one we can presently simulate, the system stops to flow 
on very long time scales and arrests in a disordered state. 

The absence of a clear diffusive regime (even at long times) and the absence of equilibrium 
conditions does not allow us to extract from the mean-squared displacement data a diffusion 
coefficient via the Einstein relation. In order to identify a relevant time scale of the structural 
dynamics, without introducing any ad-hoc model, we proceed as follows: we fix an arbitrary 
threshold value r 2 *, and consider the time r* at which r 2 (r*) = r 2 *. We repeat this analysis 
as a function of the waiting time t w . In Fig. 0] (main panel) we show our results for several 
values of r 2 *, ranging from 5 to 50 nm 2 . The data have been rescaled by t w and additionally 
shifted to maximize the overlap and stress the simple scaling with t w which, therefore, seems 
to set the only relevant time scale for diffusion. In the inset of Fig. 01 we show the raw data 
for the case r 2 * = 15 nm 2 . At short t w the data can be represented by a function of the form 
t* ~ A + Bt^, with x ~ 0.64. At longer t w we observe the linear behavior r* oc t w , expected 
for r* ~ t w . Indeed, from simple models (see, for instance, Ref. Q), it is possible to 
show that the average life-time of a potential energy landscape basin, trapping the system 
at a time of order t w , cannot exceed t w . 

Already from these results it is clear that the aging dynamics forces us to introduce a new 
relevant time scale, t w . This is of particular interest for the case of Q-dependent quantities, 
where we expect a particularly complex interplay between the length scales probed, the 
corresponding relaxation time tq and the waiting time t w . In the following Sections we 
therefore focus on the out-of-equilibrium dynamics of the density fluctuations. 
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VII. ONE-PARTICLE DYNAMICS OF DENSITY FLUCTUATIONS 



The relaxation dynamics of density fluctuations is encoded in the intermediate scattering 
function F(Q;t,t w ) = (p(Q,t)p*(Q,t w )), where the density fluctuations, p(Q,t), have been 
defined above. (We recall here that ( ) is an average over the initial conditions which amounts 
to 100 independent samples). In light scattering experiments in out-of-equilibrium condi- 
tions, information on the dynamics of the system is obtained by measuring the time autocor- 
relation function of the scattered intensity, ^(Q; t, t w ) = (I(t w + t)I(t w )). In the single scat- 
tering approximation the time autocorrelation function of the scattered intensity and the in- 
termediate scattering function are connected by the relation F(Q; t, t w ) oc ^fg^Q] t, t w ) — 1, 
where the proportionality factor depends on the detection set-up. 

We start our study of the dynamics of density fluctuations focusing on the incoherent 
(one-particle) part of the intermediate scattering function, F s (Q;t,t w ), 

1 M 

F s (Q;t,t w ) = ±(J2e~ iQ '^ i{t+tw) ~ Xi{tw) h- (10) 

i=l 

We recall here that Xi can either refer to the center of mass of platelets i, in which case 
M = N, or sites positions, in which case M = 61 x N. As already detailed above, the 
smallest value of the wave vector accessible in the simulation is Q m = 2tt/L ~ 0.055 nm" 1 . 
We have used a wave vector resolution AQ = 0.005 nm -1 and, for each value of Q, a 
maximum of 200 Q-vectors has been considered in the calculation of averages. We note 
that, although scattering experiments mainly probe the collective intermediate scattering 
function, we expect to grasp the main features of the long-time aging dynamics also from 
F S (Q; t, t w ). This choice is mainly due to the fact that one-particle dynamics is less affected 
by noise than collective dynamics and, therefore, statistics is more reliable in the former 
case. 

The quantities calculated from the position of the interaction sites also contain a dynamic 
contribution coming from rotations of the platelets; therefore, we expect the relative correla- 
tion functions to relax on time scales shorter than the analogous quantities calculated from 
the centers of mass of the platelets. Indeed, this is what we observe. In Fig. EJ we show our 
results at Q = 0.40 nm -1 , at the indicated values of the waiting time t w , both for centers of 
mass (left panel) and interaction sites (right panel). In Fig.|Hlwe show the simulation results 
at fixed t w ~ 3.5/is at the indicated values of Q. From these two figures it is evident that 
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aging is fully active in the system. In particular, the slowing down of the dynamics is larger 
the longer the aging time, and a clear Q-dependence of the relaxation is visible at constant 
t w . In the following we give a more detailed description of the effect of the aging time on 
the above dynamics. 

The most natural choice would be to fit the data to a particular model and, then, to 
study the dependence on Q and t w of the fitting parameters. In doing so we should take into 
account two important observations. First, we should note that most of the curves shown in 
Fig. El and El do not relax to zero on the available time window and this fact severely limits 
the possibility to fit the data without introducing arbitrary bounds. Second, in the present 
case t is always larger than or at most of the same order of magnitude of the waiting time 
t w . This is in contrast with the case of experiments, where it is always t w ^> t, and poses 
the question if it is reasonable to try to fit with a particular functional form a correlation 
function for times longer than the corresponding waiting time. 

The above discussion amounts to identify a reasonable relaxation time scale for the in- 
termediate scattering function which is independent of a particular model. We choose to 
use the method of the threshold we have described in the previous Section in the case of the 
mean squared displacement: we fix a particular threshold value and we consider as a relevant 
time scale the one at which the correlation function decreases from unity to the threshold 
value. Obviously, exploiting such a method we must take into account some caveats. Note, 
for instance, that the shape of the long time relaxation (possibly representable by a stretched 
exponential) clearly changes with the waiting time and, therefore, the relaxation time alone 
is not the only relevant parameter. Nevertheless, we expect to understand the correct gen- 
eral behavior and also clarify the interplay between the aging time and the length scales 
actually relevant for the aging dynamics. A-posteriori we also verify the above procedure 
in a particular case, fitting the data with a particular functional form and comparing the 
results. 

In Fig. HI we show the results for both sites (left panel) and centers of mass (right panel) 
for a threshold value / = e _1 as a function of the waiting time t w , at the indicated values of 
Q. We have checked that the value chosen for the threshold does not qualitatively change 
the results. These data confirm the picture coming from the above analysis of the mean- 
squared displacement, pointing to a strong dependence of the out-of-equilibrium relaxation 
time on the waiting time. We observe that in our data is absent the exponential dependence 
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on the waiting time of the relaxation time visible at short waiting times in experiments. 
(Note that the physical processes behind this particular feature in experiments still remains 
unexplained (l4j). 

To be more quantitative and check the analysis introduced above, we fit some of our 
data to a particular functional form, proposed some time ago and widely exploited in the 



literature (see, among others, Refs. |19|, |20j, |2l|, |42|). We choose, mainly due to a better 



quality of the data, the particular value Q = 0.7 nm 1 and fit the curves to the following 
model: 

F s (Q;t,t w ) = Ut^e-^f^ + (1 - f Q (tJ)e-W r >l < *"W lQ *' ) . (11) 

Here, we have introduced two relaxation times, i.e., T/, related to fast short-times dynamics, 
and r c , controlling the long-time relaxation. The non ergodicity parameter, fc}(t w ), and 
the stretching parameter, /3(Q;t w ), are reminiscent of the analogous quantities introduced 
for the study of the glass transition. We fit our data to Eq. (jllj) only for times tentatively 
shorter than 2t w . Fig.|H]shows the t w dependence of fitting parameters. The left panel shows 
t c , Tf and the "mean" relaxation time r m , defined as r m = t c T(1//3)//3, both for centers of 
mass of the platelets and interaction sites. All the data strongly depends on the waiting 
time and increase with t w . The right panel shows our results for (3, decreasing with t w from 
1 to about 0.4 (this behavior is also observed in experiments, see |42(), and /q, which stays 
constant at a value of about 0.1 in the case of the center of mass and increases up to 0.3 for 
sites. Fig. 121 shows the Q-dependence of the same fitting parameters at constant t w ~ 30 ns. 
Tf stays almost constant in the whole investigated Q-range, while r c shows the typical Q~ 2 



dependence [24, 



43, 



44 1 . Also, /q stays almost constant in whole considered range (data for 



platelets and sites overlap almost perfectly), while (3 shows a slight decrease with Q in both 
cases. 

We note that the above fits are particularly difficult due to the low statistics and to the fact 
that not all the curves relax to zero on the accessible time window. Moreover, the absence 
of a clear separation between time scales, makes particularly difficult the determination of 
the relaxation times. More specifically, in some cases it was also possible to fit the data with 
two comparable values for Tf and r c , obtaining a value of \ 2 comparable to the one of the 
actual fit with r/ Cr c . 
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VIII. COLLECTIVE DYNAMICS OF DENSITY FLUCTUATIONS 



The coherent intermediate scattering function at times t w and t is defined as 

1 M M 

F c (Q;t,t w ) = ^(EE^^^)' ( 12 ) 

where X{ has been defined above and S(Q) are the data included in Fig. |2l This function 
is related to the collective dynamics of the density fluctuations, and it is the quantity of 
direct relevance for the experimental measurements. In Fig. E3 we show the independence 
of F c (Q,t) at Q = 0.40 nm _1 , both for centers of mass (left panel) and sites (right panel), 
while in Fig. ^2 we show the Q-dependence of F c (Q,t) at fixed t w ~ 3.4 /is, again for the 
two cases. Also these results confirm the strong t^-dependence of the out-of-equilibrium 
dynamics and do not seem to add much more insight to the analysis of the previous Section. 
Indeed, the qualitative behavior of F c (Q;t,t w ) and F s (Q;t,t w ), as a function of both Q or 
t w , is very similar, as it is evident by comparing Fig. El to Fig. Eland Fig. ^2 to Fig. H3 
Hence, the same conclusions of Sec. IVHI can also be drawn from these results, at least at a 
qualitative level. 

In the present case we do not attempt a direct fit of the data and directly exploit the 
method of the threshold to evaluate a reasonable relaxation time scale, independent of a 
particular functional form for a fitting model. Fig. El shows our results, with a threshold 
/ = e -1 , for the data at the indicated values of Q and as a function of the waiting time, t w . 
Again, we have checked that the results do not qualitatively change with the chosen value 
of the threshold. Data are qualitatively in agreement with the results reported in Fig. [7| and 
confirm the overall picture we described in details in the previous Section. 

IX. ROTATIONAL DYNAMICS 

The strong anisotropy of the Laponite molecule implies that an important role in dynam- 
ics and structural reorganization must be played by the mutual orientations of the platelets. 
Therefore, also from the above results, we expect that the orientational behavior of the 
system, both for the single platelet and collective, must strongly couple to the dynamics of 
density fluctuations and take place on comparable time scales. This is the main motivation 
for the present Section. 
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The behaviour of the orientational degrees of freedom is accessible via the dynamics of 
the Legendre polynomials calculated from the orientation vectors of platelets. Some of those 
polynomial are directly related to the response functions detected by different experimental 
techniques. The convenient correlation functions in this case are the functions Ci(t,t w ), and 
their self-part, C\ s \t,t w ), which are defined as: 

^ N N 

Ci(t,t w ) = j^(Y,J2 Pl ^ t + t ^ ( 13a ) 

i=l 3=1 

1 N 

Cl s \t,t w ) = -(Y,Pi($i(t + t W )-tii(t«>)))- (13b) 
i=i 

Here, Ui is the normalized orientational vector of the i th molecule, pointing along the normal 
to the surface of the platelet, and Pi(x) is the Legendre polynomial of order I, with I > 1. 
Note that for I = 1 and I = 2 the functions Ci(t) can be measured in dielectric and light 
scattering experiments. We also note that it is often assumed that the cross-term in Ci(t) 
can be neglected; in this case the experiments would also yield information on Cl s \t). 

In Fig. ^1 we show our results at t w = 7 ps at the indicated values of 1 < / < 4, for the 
one-particle (left panel) and collective (right panel) cases respectively. From these data it is 
evident that the smaller the value of I, the longer the time scale of the relaxation dynamics. 
In particular, note that in the collective case, the curves seem to saturate to a finite value 
at long time, pointing to a progressive freezing of the orientational degrees of freedom at 
very long times. Unfortunately the limited total time extension of the present computer 
simulation does not allow us to be conclusive on this point. 

In Fig. we plot our results for I = 2 (panels a) and b)) and 4 (panel c) and d)) at the 
indicated waiting times, both for the self (left panels) and collective (right panels) cases. 
Again, also in the case of the rotational dynamics, a strong aging time dependence is evident 
and at longer waiting times the rotational dynamics slows down significantly. 

A threshold analysis similar to the one described above for the dynamics of density 
fluctuations allows one to identify a relaxation time. Our results are shown for 1 < I < 4 
in Fig. where the threshold has been fixed to a value / = 0.8. The overall qualitative 
behavior is analogous to the one found in the case of the translational dynamics, and the time 
scales, except for the case 1 = 1, are also comparable in the two cases. Concluding, our data 
support a strong coupling between translational and orientational degrees of freedom 46^ . 
as expected for a strongly anisotropic molecule. 
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X. CONCLUSIONS 



In this work we have presented a systematic study of the out-of-equilibrium and aging 
dynamics of a model for a Laponite colloidal suspension. The choice of the system considered 
has been mainly suggested by the terrific amount of experimental work already present 
in the literature. A realistic site-site, purely repulsive Yukawa-like interaction potential 
between Laponite platelets, has been considered. We have used the Brownian dynamics 
technique for extensive molecular dynamics simulations. In particular, we have proposed a 
new algorithm which, although neglecting hydrodynamic effects, correctly takes into account 
the strong anisotropy of the Laponite platelet, introducing two different friction coefficients, 
respectively parallel and perpendicular to the direction of the symmetry axis of the platelet. 

An extensive ensemble of independent systems has been driven out-of-equilibrium by 
an instantaneous jump from a high temperature to room temperature at constant volume 
fraction and realistic configurations have been stored for the analysis. In particular, we 
have focused on the long-time aging dynamics, calculating the mean-squared displacement 
and the intermediate scattering functions. Both coherent and incoherent dynamics have 
been studied in details. The dynamics of the density fluctuations has been found to be 
strongly dependent on the aging time and a qualitative discussion has been performed on 
the interplay between the different time and length scales which play a significant role in 
the relaxation. A study of the rotational dynamics has completed the study of the dynamics 
of the Laponite colloidal suspension. The main result is that aging dynamics also strongly 
affects the orientational degrees of freedom, which relax on time scales comparable to the 
ones typical of translational modes. 

The study presented here strongly suggests that, indeed, in the absence of attractive 
interactions between molecules, a disordered arrested state can be generated by the only 
presence of Yukawa repulsive interactions. The long range of the repulsive interaction makes 
it possible to arrest the dynamics of the colloidal suspension even if the effective packing 
fraction is only a few percent. In this respect, cages are defined not by the physical size of 
the hard-core but by the effective range of the screened electrostatic potential. While in the 
case of spherical particles, Yukawa interactions alone hardly generate a glass, due to the fast 
rates of crystallization induced by the long range of the interaction, the anisotropy of the 
platelets appears to be able to self-generate a sufficient local disorder which stabilizes the 
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metastability of the arrested disordered state. 

In this respect, despite our data are not conclusive, due to the difficulty of extending 
the time length of the simulation, charged platelets may give rise to Wigner glasses. Future 
work, in line with the present approach, must focus on the role played in the dynamics by 
both the salt concentration and the presence of opposite sign rim charges. 
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FIGURE CAPTIONS 



FIG 1: Relaxation of potential energy as a function of time after the temperature jump. 
100 independent realizations of the system have been equilibrated at high temperature, 
T = 500 K, and, at t w = 0, instantaneously quenched to T = 300 K. In the plot we 
show the data for each run (lines) and the averaged value (symbols). Note that the energy 
stays positive, due to the purely repulsive nature of the interaction site-site potential. The 
relaxation of the energy follows a power law behavior of exponent 7 ~ 0.35 and the value 
predicted at equilibrium (estimated at Eoo ~ 206.6 eV) is not reached on the total simulation 
time scale. Therefore, the systems stay out-of-equilibrium for the total duration of the 
present simulation. 

FIG 2: Main Panel: Equilibrium static structure factor, S(Q), calculated both for the 
platelets centers of mass (open circles) and interaction sites (closed circles). Each point is 
an average over 100 system configurations for t > 10 6 ps. The smallest value of Q accessible 
in our simulation is Q m = 2n/L ~ 0.05 nm -1 . For the wave- vector average, we have chosen 
a resolution AQ = 0.005 nm -1 and up to 200 Q-vectors have been considered for each value 
of Q. Error bars are smaller than the symbols. Inset: Variation of the static structure 
factor calculated on sites at short times at the indicated values of t. A slight change in the 
structure of the second sub-peak is evident. 

FIG 3: Mean squared displacement (r 2 (t,t w )), calculated for the platelets centers of mass, 
as discussed in the text, as a function of time for different waiting times t w . Shorter waiting 
times are on the top. The initial transient dynamics at short t is followed by a plateau whose 
life-time increases on increasing t w . Next the curves move from the plateau, saturating to a 
(nearly) constant value at very long times. From these data we could argue that the system 
eventually stops to flow on very long time scales. 

FIG 4: Relevant dynamic time scale r* calculated from the mean squared displacement, 
as discussed in the text. Main panel: r* as a function of t w calculated for different values 
of r 2 *, ranging from 5 to 50 nm 2 . Data have been rescaled by t w and additionally shifted 
to maximize the overlap and stress the simple scaling with t w , which sets the relevant time 
scale. Inset: r* as a function of t w for r 2 * = 15 nm 2 . At short t w data can be represented by 
a power law behavior, r* ~ A + Bt^, with x ~ 0.64 (solid line); at longer t W) one observes 
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the expected linear behavior r* oc t w (dashed line) 
FIG 5: t^-dependence of the incoherent (one-particle) intermediate scattering functions 
F s (Q;t,t w ) at Q = 0.40 nm" 1 , at the indicated values of t w . a) Data calculated from the 
positions of the platelets centers of mass, b) Data calculated from sites, also containing 
orientational dynamics contributions. The effect of aging is evident from these data. 
FIG 6: Q-dependence of the one-particle intermediate scattering function F s (Q;t,t w ) at 
t w = 3448933 ps at the indicated values of Q. a) Data calculated on the platelets centers 
of mass, b ) Data calculated on interaction sites. A strong dependence on Q is evident from 
these data. 

FIG 7: Relaxation time for the incoherent intermediate scattering functions calculated by 
the threshold method with / = e _1 for both centers of mass and sites as a function of the 
waiting time t w . The indicated values of Q have been considered. 

FIG 8: Waiting time dependence of the parameters of the fit discussed in the text for the 
incoherent scattering function at Q = 0.7 nm -1 . 

FIG 9: Parameters of the fits discussed in the text for the incoherent scattering function 
at t w = 29390 ps as a function of the wave vector Q. 

FIG 10: ^-dependence of the coherent (collective) intermediate scattering function, 
F c (Q;t,t w ), at Q = 0.40 nm -1 , at the indicated values of t w (longer t w 's on top). Data 
calculated on the platelets centers of mass (a)) and on interaction sites (b)) are qualitatively 
in agreement. 

FIG 11: Q— dependence of the coherent intermediate scattering function, F c (Q;t,t w ), at 
t w = 3448933 ps at the indicated values of Q. Data calculated on the platelets centers of 
mass (a)) and on interaction sites (b)) show a qualitative similar behavior. 
FIG 12: Relaxation time for the coherent intermediate scattering functions calculated by 
the threshold method with / = e -1 for centers of mass and sites as a function of the waiting 
time t w at the indicated values of Q. 

FIG 13: Correlation functions of the Legendre polynomials of the order 1 < I < 4 at 
t w ~ 7 ps. We show the self contribution (left panel) and collective total case (right panel). 
The orientational relaxation dynamics is slower the higher the value of the order I. 
FIG 14: Top: Waiting time dependence of the orientational correlation functions for I = 2 
(longer waiting times on the top), both for the self (left panel) and collective (right panel) 
cases. Bottom: Waiting time dependence of the orientational correlation functions for I = 4 
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(shorter waiting times on the top) both for the self (left panel) and collective total case 
(right panel). 

FIG 15: Aging time dependence of the orientational relaxation time for the correlation 
functions of the Legendre polynomials of order 1 < I < 4 determined by means of the 
threshold method, as discussed in the text. The threshold has been fixed to the value of 
/ = 0.8. Left panel: Data calculated from the self part of the correlation functions; Right 
panel: Data calculated from the total collective correlation functions. 



23 




24 




25 



!""' 1 










— 397 




— 1678 




— 7029 




— 29374 = 




— 122721 : 




— 512646 , 




— 3448930 



' 1 2 1 4 S (S 7 

10 10 10 10 10 10 10 

1 - 1 „ [PS] 



FIG. 3: 



26 




io 1 io 2 io 3 io 4 io 5 io 6 io 7 
FIG. 4: 



27 




FIG. 5: 



28 




io 2 io 3 io" io 5 io 6 io 7 io 3 io 4 io 5 10* io 7 
t-t[ps] 



FIG. 6: 



29 




10 1 io 2 io 3 io 4 io 5 io" 10 1 io 2 10 3 1() 4 10 5 10 6 



FIG. 7: 



30 



10 i o 1 v cms 1 "77 

I □ T c , CMs ^ 

10 7 i ♦ t m , Sites n - 

: . x c , Sites \ n 

' • T p Sites 

10 S = o t CMs .oS8 a 



4 Q = 0.7[nm- 1 ] 



ggSooooooo 



1 111 

8 


' 

b) 


8 








8 .°o 


- 


'•So 

. p, Sites 





" o p, CMs 




- . /q, Sites 


oo 


- □ / Q ,CMs 




IHHBHHlHB'nn 


□ □ □ aanP - 







10° 10* 



10 



FIG. 8: 



31 




10"' 10° 10"' 10 1 10° 

Q [nm '] 



FIG. 9: 



32 




FIG. 10: 



33 




FIG. 11: 



34 




FIG. 12: 



35 




io 2 io 3 io" io 5 io 6 io 7 io 3 io 4 io 5 10* io 7 
t-t w [ps] 

FIG. 13: 



36 




a) Self b) Collective 



10 10 10 10 
t-t [ps] 



10" 



10' 




io 2 io 3 io 4 io 5 io 6 io 7 io 3 io 4 io 5 10* io 7 
t-t[ps] 



FIG. 14: 



37 



1 1 1 1 1 • 


»"l 1 1 1 1 1 1 

_/// ; 












b) Collective 





io 1 io 2 io 3 io 4 io 5 io' io 2 io 3 io 4 io 5 io 6 io 7 



FIG. 15: 



38 



